#ifndef FLUX_
#define FLUX_

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <vector>
#include <map>
#include <cmath>
#include "ANN/ANN.h"
#include <omp.h>

#include "DGMatrix.h"
#include "DGIntMatrix.h"
#include "mesh.h"
#include "saveGmsh.h"
#include "initialSol.h"
#include "functionSpace.h"
#include "muparser/muParser.h"

using namespace std;

enum PType    {EQP,ENP,FQP,FNP,BQP,BNP};
enum BCType   {Periodic,Dirichlet};
//------------------------------------------------------------------------------
class Flux{
  
private:
  
  void getBondaryConditions(ClassParameter&);
  
protected:
  
  int             nbProcs;
  int             nbEq;
  Mesh*           mesh;
  FsContainer*    fs;
  BCType          BCtype;
  DGMatrix        UBNP;                                                         //(nbFNP,nbB)

public:
  
  Flux(Mesh*,FsContainer*,ClassParameter&);
  ~Flux(){};
  
public:
  
  virtual void getFlux           (int,const DGMatrix*,DGMatrix*,PType,int,int)=0;
  virtual void getFluxDerivative (int,const DGMatrix*,DGMatrix*,PType,int,int)=0;
  virtual void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                                  const DGMatrix*,const SolContainer*,const SolContainer*,
                                  DGMatrix*,const DGVector*,int)=0;

  virtual double getBetaMax()=0;
  virtual double getBetaMin()=0;
  virtual double getMaxBetaOverHmax()=0;
  virtual double getMaxBetaOverHmin()=0;
  
  virtual double getBCValue(int bID,int nI){return UBNP(nI,bID);}
  
};
//------------------------------------------------------------------------------
class LinearTransport:public Flux{
  
protected:

  SolContainer  BetaEQP;
  SolContainer  BetaFQP;
  SolContainer  BetaENP;
  SolContainer  BetaFNP;
  SolContainer  BetaBQP;
  SolContainer  BetaBNP;
  
  double BetaMax;
  double BetaMin;
  double maxBetaOverHmax;
  double maxBetaOverHmin;
  
public:
  
  LinearTransport(Mesh* mymesh,FsContainer* myfs,ClassParameter& param);//:Flux(mymesh,myfs,param);
  ~LinearTransport(){};
  
public:
  
  void         getFlux           (int,const DGMatrix*,DGMatrix*,PType,int,int);
  void         getFluxDerivative (int,const DGMatrix*,DGMatrix*,PType,int,int);
  virtual void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                                  const DGMatrix*,const SolContainer*,const SolContainer*,
                                  DGMatrix*,const DGVector*,int)=0;  
  
  double getBetaMax()         {return BetaMax;}
  double getBetaMin()         {return BetaMin;}
  double getMaxBetaOverHmax() {return maxBetaOverHmax;}
  double getMaxBetaOverHmin() {return maxBetaOverHmin;}
  
};
//------------------------------------------------------------------------------
class UpwindLinearTransport:public LinearTransport{
  
public:
  
  UpwindLinearTransport(Mesh* mymesh,FsContainer* myfs,ClassParameter& param):LinearTransport(mymesh,myfs,param){};
  ~UpwindLinearTransport(){};

public:
  
  void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                          const DGMatrix*,const SolContainer*,const SolContainer*,
                          DGMatrix*,const DGVector*,int);  

};
//------------------------------------------------------------------------------
class CenterLinearTransport:public LinearTransport{
  
public:
  
  CenterLinearTransport(Mesh* mymesh,FsContainer* myfs,ClassParameter& param):LinearTransport(mymesh,myfs,param){};
  ~CenterLinearTransport(){};
  
public:

  void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                          const DGMatrix*,const SolContainer*,const SolContainer*,
                          DGMatrix*,const DGVector*,int);  
  
};
//------------------------------------------------------------------------------
class KPP:public Flux{
  
protected:
  
public:
  
  KPP(Mesh* mymesh,FsContainer* myfs,ClassParameter& param);
  ~KPP(){};
  
public:
  
  void         getFlux           (int,const DGMatrix*,DGMatrix*,PType,int,int);
  void         getFluxDerivative (int,const DGMatrix*,DGMatrix*,PType,int,int);
  virtual void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                                  const DGMatrix*,const SolContainer*,const SolContainer*,
                                  DGMatrix*,const DGVector*,int)=0;  
    
  double getBetaMax()         {return 0;}
  double getBetaMin()         {return 0;}
  double getMaxBetaOverHmax() {return 0;}
  double getMaxBetaOverHmin() {return 0;}
  
};
//------------------------------------------------------------------------------
class UpwindKPP:public KPP{
  
public:
  
  UpwindKPP(Mesh* mymesh,FsContainer* myfs,ClassParameter& param):KPP(mymesh,myfs,param){};
  ~UpwindKPP(){};
  
public:
    
  void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                          const DGMatrix*,const SolContainer*,const SolContainer*,
                          DGMatrix*,const DGVector*,int);  
  
};
//------------------------------------------------------------------------------
class CenterKPP:public KPP{
  
public:
  
  CenterKPP(Mesh* mymesh,FsContainer* myfs,ClassParameter& param):KPP(mymesh,myfs,param){};
  ~CenterKPP(){};
  
public:
    
  void getNumericalFlux  (const DGMatrix*,const SolContainer*,const SolContainer*,
                          const DGMatrix*,const SolContainer*,const SolContainer*,
                          DGMatrix*,const DGVector*,int);  
  
};
//------------------------------------------------------------------------------
class FluxFactory{
public:
  static Flux* CreateNewFlux(Mesh*,FsContainer*,ClassParameter&);
};
#endif
